function [u] = ex32_SolveSystem_sol(dp_lme,s_near,x_nodes,x_samples,w_samples,parameters,options)
% function [u] =
% ex32_SolveSystem_sol(dp_lme,s_near,x_nodes,x_samples,w_samples,parameters,options)
%
% This function assembled the stiffness matrix K and the right hand side
%     rhs corresponding to the cantilever beam problem explained in
%     Timoshenko's book.
%
% Input:
%    p_lme     : shape functions
%    s_lme     : shape functions gradient
%    s_near   : list of neighbors
%    x_nodes   : node points
%    x_samples : sample point
%    w_samples : gauss weigth for each sample point
%    parameters: L (length), D (diameter), nu (Poisson coefficient), E
%                (Young modulus)
%    options   : lme options
%
% Output:
%    u     : vectorial displacement field
%


% Material parameters
E  = parameters.E;
nu = parameters.nu;

nPts = size(x_nodes,1);
sPts = size(x_samples,1);


%% ------------------------------------------------------------------------
% The right hand side rhs is computed
[rhs ind_Dirichlet] = ex32_RHS_sol(x_nodes,options,parameters);


%% ------------------------------------------------------------------------
%  The stiffness matrix is assembled 
nn=0;
for k=1:sPts
  nn = max(nn, length(s_near{k}));
end
nn = min(nn,nPts);

%K = spalloc(2*nPts,2*nPts,2*nn*nPts);
K = zeros (2*nPts,2*nPts);
C_stiff=E/(1+nu)/(1-2*nu)*[1-nu,   nu,          0 ;...
	                           nu, 1-nu,          0 ;...
	                            0,    0, (1-2*nu)/2];

                            
                            
%% Assembly K matrix per element - constitutive part

   % Ortiz 2010, "Maximum-entropy meshfree method for compressible and
   % near-incompressible elasticity"
   
   % Ensamblado ecuación 26, matriz de rigidez modificado
   

for ig=1:sPts

   nact=length(s_near{ig});
   active=s_near{ig};
   
   for i=1:nact
       
   act = active(i);
   
   % Calculo la matriz Ba correspondiente al nodo (act) a en el punto (ig).
   % Ecuación (24d)
   
   B(1,1,act) = dp_lme{ig}(i,1);
   B(1,2,act) = 0;
   B(2,1,act) = 0;
   B(2,2,act) = dp_lme{ig}(i,2);
   B(3,1,act) = dp_lme{ig}(i,2);
   B(3,2,act) = dp_lme{ig}(i,1);
   
   % Nodos cercanos al punto de Gauss ig
   
   %elem_nodes = samplesAdjacency_fem (x_gauss,conectivity,gPts);
   
   % Elementos para modificación del tensor B
   
   % Matriz Na
   
   %
   
   % Puntos de Gauss sobre los bordes
   
   
   
   % 

   % Calculo la matriz Ga
   
   % 
   
      
     
   end

  %assembly
  
  for i=1:nact
      for j=1:nact
          i_act = active(i);
          j_act = active(j);
          ix = 2*i_act-1;
          iy = 2*i_act;
          jx = 2*j_act-1;
          jy = 2*j_act;
          
          K_ig = B(:,:,i_act)' * C_stiff * B (:,:,j_act);
          
          K (ix,jx) = K (ix,jx) + K_ig(1,1) *w_samples(ig);
          K (iy,jy) = K (iy,jy) + K_ig(2,2) *w_samples(ig);
          K (ix,jy) = K (ix,jy) + K_ig(1,2) *w_samples(ig);
          K (jy,ix) = K (ix,jy);
          
      end
  end
  
 
end

%% Assembly K matrix per element - Pressure part






%% Dirichlet BCs are applied

%(x=0,y=0)  ux=0 uy=0
ind = ind_Dirichlet(1);
K(2*ind,:)         = 0;
K(:,2*ind)         = 0;
rhs(2*ind)         = 0;
K(2*ind,2*ind)     = 1;
K(2*ind-1,:)       = 0;
K(:,2*ind-1)       = 0;
rhs(2*ind-1)       = 0;
K(2*ind-1,2*ind-1) = 1;

%(x=0,y=D/2)  ux=0
ind = ind_Dirichlet(2);
K(2*ind-1,:)       = 0;
K(:,2*ind-1)       = 0;
rhs(2*ind-1)       = 0;
K(2*ind-1,2*ind-1) = 1;



%% ------------------------------------------------------------------------
% The system is solved
u=K\rhs;


%% Assembly original

% for ig=1:sPts
% 
%    nact=length(s_near{ig});
%    active=s_near{ig};
%    
% 
%    
%  
%   B_ig=zeros(3,2*nact);
%   B_ig(1,1:2:2*nact)=dp_lme{ig}(:,1)';
%   B_ig(2,2:2:2*nact)=dp_lme{ig}(:,2)';
%   B_ig(3,2:2:2*nact)=dp_lme{ig}(:,1)';
%   B_ig(3,1:2:2*nact)=dp_lme{ig}(:,2)';
%   K_ig_loc=B_ig'*C_stiff*B_ig;
%   
%   
%   %assembly
%            
%           
%   K(2*active(:)-1,2*active(:)-1) = ...
%       K(2*active(:)-1,2*active(:)-1) + ...
%       K_ig_loc(1:2:2*nact,1:2:2*nact)*w_samples(ig);
%   K(2*active(:),2*active(:)-1) = ...
%       K(2*active(:),2*active(:)-1) + ...
%       K_ig_loc(2:2:2*nact,1:2:2*nact)*w_samples(ig);
%   K(2*active(:)-1,2*active(:)) = ...
%       K(2*active(:)-1,2*active(:)) + ...
%       K_ig_loc(1:2:2*nact,2:2:2*nact)*w_samples(ig);
%   K(2*active(:),2*active(:)) = ...
%       K(2*active(:),2*active(:)) + ...
%       K_ig_loc(2:2:2*nact,2:2:2*nact)*w_samples(ig);
% end

